pwd
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar
import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
from scipy.stats import boxcox, probplot
bios = gpd.read_file('biomes_continents')
bios.columns
area = bios['SUM_Area']*1e-6
data_columns =bios.columns[4:4+17*12]
inter_cols = {}
for yr in range(17):
inter_cols[yr] = data_columns[12*yr:12*yr +12]
intra_cols = {}
for mo in range(12):
intra_cols[mo] = [data_columns[mo + yr *12 ] for yr in range(17)]
for yr in range(17):
inter = bios[ inter_cols[yr]].sum(axis=1)/area
label = 'inter_' + str(yr +2001)
bios= bios.assign(label=inter)
bios = bios.rename(columns={'label': label})
inter_columns = bios.columns[-17:]
inter_max = bios[inter_columns].max().max()
for mo in range(12):
intra = bios[ intra_cols[mo]].sum(axis=1)/(17*area)
label = 'intra_{:02}'.format(mo)
bios = bios.assign(label=intra)
bios = bios.rename(columns={'label': label})
intra_columns = bios.columns[-12:]
intra_max = bios[intra_columns].max().max()
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100):
fig, ax = plt.subplots(figsize=(12, 8))
ax.set_aspect('equal')
ax = df.plot(ax =ax, column=column, cmap='OrRd',vmin=vmin, vmax=vmax)
ax.set_axis_off()
norm = cls.Normalize(vmin, vmax)
cmmapable = cm.ScalarMappable(norm, 'OrRd')
cmmapable._A = []
cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal")
cb.set_label(title, fontsize=20, family='Times New Roman')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')
plt.show()
fig.savefig(name, dpi=dpi, bbox_inches='tight')
lambdas = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
lambdas
probplot(bios.inter_2017, plot=plt.gca());
probplot(test_column, plot=plt.gca());
maxs = []
mins = []
for yr in range(17):
column= 'inter_{}'.format(2001 + yr)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
for yr in range(17):
column = 'inter_{}'.format(yr +2001)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios = bios.assign(test_column =test_column)
yr_name = str(yr +2001)
#"Burnt area pixels (#/km2) by biome in 2001
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +yr_name
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'inter_biomes_parts_index_{}'.format(yr +2001)
draw_map(bios,'test_column',title=title,name=name)
lambdas = []
for mo in range(12):
column= 'intra_{:02}'.format(mo)
test_column, lbd = boxcox(bios[column] + 1e-25)
lambdas.append(lbd)
lmbda= sum(lambdas)/len(lambdas)
lambdas
probplot(bios.intra_11, plot=plt.gca());
probplot(test_column, plot=plt.gca());
maxs = []
mins = []
for mo in range(12):
column= column= 'intra_{:02}'.format(mo)
test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
maxs.append(test_column.max())
mins.append(test_column.min())
box_min = min(mins)
box_delta =max(maxs) - box_min
for mo in range(12):
column = 'intra_{:02}'.format(mo)
test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
bios= bios.assign(test_column =test_column)
mo_name = calendar.month_name[mo + 1]
#"Burnt area pixels (#/km2) by biome in Janury (2001-2017)"
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +mo_name + ' (2001-2017)'
#title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by climate 2001-2017 ' + mo_name
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'intra_biomes_parts_index_{:02}'.format(mo)
draw_map(bios,'test_column',title=title,name=name)
total = bios[data_columns].sum(axis=1)/17
total = total/area
bios= bios.assign(total=total)
column = 'total'
test_column, lbd = boxcox(bios[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
bios = bios.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in (2001-2017)'
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes_parts_index'
draw_map(bios,'test_column',title=title,name=name)